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ABSTRACT 

We investigate the dynamical interaction of a central star cluster surrounding a super-massive 
black hole and a central accretion disk. The dissipative force acting on stars in the disk leads 
to an enhanced mass flow towards the super-massive black hole and to an asymmetry in the 
phase space distribution due to the rotating accretion disk. The accretion disk is considered as 
a stationary Keplerian rotating disk, which is vertically extended in order to employ a fully self- 
consistent treatment of stellar dynamics including the dissipative force originating from star-gas 
ram pressure effects. The stellar system is treated with a direct high-accuracy TV-body integration 
code. A star-by-star representation, desirable in TV-body simulations, cannot be extended to real 
particle numbers yet. Hence, we carefully discuss the scaling behavior of our model with regard 
to particle number and tidal accretion radius. The main idea is to find a family of models for 
which the ratio of two-body relaxation time and dissipation time (for kinetic energy of stellar 
orbits) is constant, which then allows us to extrapolate our results to real parameters of galactic 
nuclei. Our model is derived from basic physical principles and as such it provides insight into the 
role of physical processes in galactic nuclei, but it should be regarded as a first step towards more 
realistic and more comprehensive simulations. Nevertheless, the following conclusions appear 
to be robust: the star accretion rate onto the accretion disk and subsequently onto the super- 
massive black hole is enhanced by a significant factor compared to purely stellar dynamical 
systems neglecting the disk. This process leads to enhanced fueling of central disks in active 
galactic nuclei and to an enhanced rate of tidal stellar disruptions. Such disruptions may produce 
electromagnetic counterparts in form of observable X-ray flares. Our models improve predictions 
for their rates in quiescent galactic nuclei. We do not yet model direct stellar collisions in the 
gravitational potential well of the black hole, which could further enhance the growth rate of the 
black hole. Our models are relevant for quiescent galactic nuclei, because all our mass accretion 
rates would give rise to luminosities much smaller than the Eddington luminosity. To reach 
Eddington luminosities, outflows and feedback as in the most active QSO's other scenarios are 
needed, such as gas accretion after galaxy mergers. However, for AGNs close to the Eddington 
limit this process may not serve as the dominant accretion process due to the long timescale. 

Subject headings: accretion disks - galaxies, galactic nuclei - celestial mechanics - stellar dynamics, N- 
body simulations 
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1. Introduction 

Kinematic and photometric data of galactic nu- 
clei have revealed that super-massive central black 
holes (S MBHs) are ubiquitous in most galaxies 
(see e.g. iKormendv fe Richstonej 19951 for a re- 
view). Detailed information on their photometric 
profiles and shapes of spectral line profiles allow 
us, in certain limits, to deduce the true shape of 
the distribution function in phase space of such 
systems. The distribution function determines 
the rate at which stars will come close to the 
SMBH, be tidally disrupted or destroyed by direct 
collisi ons and eventually accreted onto the blac k 
hole dFrank fc Reedll976l : iBahcall fe Wohlll976l) . 
Quasars, however, which are the most luminous 
witnesses of accretion activity onto SMBHs are 
observed already in the young universe at red- 
shifts of z > 5. It is difficult to explain how black 
holes can grow so quickly to the observed high 
masses (10 6 — 1O 9 M ) by pure star accretion. It 
has been therefore argued that massive seed black 
holes form already by dissipative and viscous col- 
lapse, possibly accompanied by the formation of 
massive stars and their coalescence, at the time of 



galaxy formation (lColgatdll 967; Spitzcr & Saslaw 



1966[ ISpitzer fc Stond Il967t ISandersI Il970i: iRees 
1984)). In the hierarchical picture of galaxy for- 
mation, the most massive dark halos with small 
angular momentum can account for the early for- 
mation of the most massive black holes. The 
SMBH forms by a centrally focused collapse en- 
tering a phase of a dense super-massive gaseous 
object which is supported by star-gas i nteractions 
dBisnovatvi-Kogan fc Svunvaevl Il972t [Vilkoviskiil 
19751 lHaralll978t lLangbein et al.lll990l ). The in- 
teraction of a compact stellar cluster with a mas- 
sive central object in the fo rm of a supersta r 
or SMBH was con s idered by IVilkoviskiil (|l975h : 
Hilisl (|l975h : lHaral (jl978l ). The evolution of the 
dense non-rotating stel l ar cluster was studied b; 
Spitzer fc Saslawl (| 1 9661 ) ; iBisnovatvi-Koganl (j 1 9 7 
among others, and the evolution o f the g as sphere 
was considered bv lLangbein et all (|l990h . 

These studies have commonly neglected angular 
momentum, which should not be neglected during 
mergers nor in the intrinsic structure of galaxies. 
Spectra of active galactic nuclei suggest that there 
should be gas in the form of a massive central ac- 
cretion disk (AD) in which all interstellar matter 



settles before it may be accreted. The origin of 
the AD could be inflowing cool gas during mergers 
and/or debris by direct stellar collisions or stellar 
evolutionary processes, depen ding on the evolu- 



tionar y state of the host galaxy. lArtvmowicz et al 



( 19931 ) provide a theoretical model framework in 
which star-gas interactions and the build-up of 
massive stars generate a massive AD around the 
central SMBH. 

Recent work shows that gaseous disks in galac- 
tic centers around SMBHs are important for the 
dynamics and the morphology of central regions of 
galaxies and for the evolution of single or multi- 
ple central black hol e s in many ways. For example 
ICuadra et all (<2009l ): ICallegari et al.1 (<2009l . l201lh 
study the role of small scale disks for the accelera- 
tion of binary black hole mergers, through dissipa- 
tion of k inetic energy in the di sk, after a galactic 
merger. iBaruteau et al. ( 2011 ) discuss the hard- 
ening of stellar binaries in circumnuclear disks and 
their subsequent interactions with central black 
holes, which may lead to high velocity stellar es- 
capers. 



Shakura fc Sunvaevi (|1973l ) developed a model 



for a stationary AD, which has been the basis of 
many investigations since then. Close to the inner 
boundary of the AD, at a few Schwarzschild radii 
rs, gener al relativistic effects must be taken into 
account. iNovikov fc Thornd (jl973l ) extended the 
standard disc model in that regime. 

If the inner part of the AD reaches a crit- 
ical surface density the interaction of orbits 
with the gas (or more correctly energy and 
momentum transfer of stars due to ram pres- 
sure effects, henceforth denoted as dissipation) 
cannot be neglected anymore. Stellar interac- 
tions with accretion disks were also considered 



by IV ilkoviski i fc Bekbosarovl (|1982I ): IVilkoviskij 
(|l983l ): ISver et al.l (|l99ll ). More detailed inves 



tigations of the stellar or bits, crossing accretion 
disks were presented in I Vokrouhlickv fc Karas 



1998) and i nlSubr fc Karad(ll999l) . ISubr fc Ka~ 



19991) and ISubr et al.l (120041) assumed an in- 
finitely thin disk interacting with stars; they found 
that the interaction between the disk and the 
stars (star-gas dissipation) will deplete counter- 
rotating stars, create a flattening of the large-scale 
structure of the system and initiate anisotropies 
(i.e. changes of the eccentricity distributions). 
These studies did not take into account any feed- 
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back onto the disk or any finite-thickness effects 
of the disk. They tried to find a stationary state 
in which the transport of stars into the central 
galactic region is balanced by removal of stars 
from the central disk. Evolution timescales and 
initial conditions to reach this equilibrium was 
missing. 

Around SMBHs the AD may extend to the 
parsec scale, orders of magnitude larger than r§ . 
Rauch (I1995L Il999ft studied in detail the impact 
of the AD on the orbits and the distribution of 
stars with test particle simulations. Their work 
includes relativistic effects, even for the case of ro- 
tating Kerr black holes. They find that the orbital 
inclination with respect to the AD declines quickly 
as soon as the dissipative force becomes effective. 
Additionally they find a steepening of the central 
density cusp, due to the combined effect of rela- 
tivistic orbit migration and stellar collisions. 

Th e semi-analytic model bv lVilkoviskii fc Czernv 
(|2002| ) raises the point that two-body relaxation of 
the stars within and near the disk will tend to ele- 
vate trapped stars again out of the disk, and that 
the competition between relaxation and dissipa- 
tion will define a stationary state of the system, 
with some well-defined stati onary flux of stars go- 
i ng do wn to the black hole. IVilkoviskii fc Czernv 
(|2002l ) compared the star-star two-body interac- 
tions with the star-disk interactions and concluded 
that the latter is stronger in the inner parts of the 
accretion disk and vice versa in the outer parts 
based on nearly circular orbits of the stars. They 
derived analytic approximations for the effective 
inclination of stellar orbits, where the inflow to the 
SMBH takes place. They also derived a critical 
radius inside which direct stellar collisions must 
be taken into account. 

All of the mentioned papers so far have both 
strong and weak points. iRauchl d!999h is the only 
paper to combine general relativity effects and 
stellar collisions, but considered only the central 
region where the (spherically symmetric) poten- 
tial is dominated by the super-massive black hole. 
They used an approximate model of stellar dy- 
namics based on a Monte Carlo technique, which 
requires a sphe r ically symmetric central star clus- 
ter. In Rauch ( 19951 ) they combined relativistic 
effects with star-disk interactions, but the disk 
is assumed to be infinitely thin. The infinitely 
thin disk approximation was also used in other 



work bvlSubr fe Karasl (fl999h : ISubr etaD t004): 
Vilkoviskii fc Czernvl (|20olV ). 



Fol lowing the ideas of IVilkoviskii fc Czernv 
(|2002f ) we investigate by self-consistent direct TV- 
body simulations, the interaction of a central com- 
pact stellar cluster (CSC) with the AD in active 
galactic nuclei. We focus on the mutual interplay 
of two-body relaxation and the depletion of stars 
by the dissipative force in the AD as a secular 
long-term evolution of the stellar mass distribution 
and the velocity distribution function of the CSC. 
In this paper we report results on a new model of 
star-disk interactions in galactic nuclei. Our focus 
lies on the correct and accurate representation of 
the stellar orbital motion crossing the disk, by im- 
plementing a disk with its density distribution in 
a full three-dimensional iV-body simulation. We 
have added the force and force time derivative 
to the standard Hermite scheme (see below) as a 
function of local disk density and velocity in three 
dimensions. This is the most general approach, 
and later it will allow us to include evolving mod- 
els of the AD and appropriately model the mass 
and energy transfer between the AD and CSC. 

Our first results concern the enhanced accre- 
tion rate onto the central SMBH due to the in- 
teraction with the AD, in the regime where the 
relaxation timescale is comparable to the dissipa- 
tion timescale. This is, to our knowledge, the first 
approach to study the competition between relax- 
ation and star-gas dissipation in a direct simula- 
tion. Since the direct simulations are not yet able 
to reach realistic particle numbers and spatial res- 
olution, we will perform a careful scaling analysis 
to show in which way our numerical simulations 
have to be interpreted for real astrophysical galac- 
tic nuclei. 

This paper is organized as follows: in Sect. [2] 
we describe the accretion disc and the dissipative 
force in detail, Sect. |3] gives the numerical real- 
izations of the system, in Sect. 0] the results are 
discussed and in Sect. [5] a summary and conclu- 
sions are presented. 

In follow-up papers we will discuss the depen- 
dence of the dissipation on the orbital param- 
eters and the phase space evolution of the cusp 
stars in detail and take into account the feedback 
of the star-gas interaction on the AD properties. 
Detailed studies of migration of stars, binaries, 
and black holes inside the disk towards the cen- 
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ter, and its observational consequences will also 
be included in future work. 

2. Physics of the accretion disk and the 
star-gas interactions 

We consider an AGN model consisting of three 
subsystems: i) a compact stellar cluster (CSC) 
with mass M c \ describing the inner part of the 
galactic center. It is spherically symmetric, non- 
rotating, and in dynamical equilibrium, ii) An 
accretion disk (AD) with mass Md = HdM c \. It is 
vertically extended, non-evolving, and has a Ke- 
plerian rotation curve, iii) A central super mas- 
sive black hole (SMBH) with mass M^h = Mbh-^ci- 
The motion of each star m\ of the CSC is deter- 
mined by the mutual gravitational interaction of 
the stars, the gravitational force of the SMBH, 
and by a dissipative force Fa from the AD. The 
equation of motion is given by 



where R 2 



Gr 



GM hh ?i , F d 



(1) 



where fy = fj — f] with n, rj the positions of stars i 
and j, respectively. Since the AD has a small mass 
compared to the black hole and the enclosed stellar 
cluster, we neglect the gravitating effect of the disk 
on the system. Numerical details for computing 
the forces are given in Sect. [3] 

2.1. The accretion disk 

We are interested in the dynamical action of 
the AD on the stellar component. Therefore we 
adopt a 3-dimensional axisymmetric stationary 
disk model, which is differentially rotating with 
the local circular velocity. For the inner struc- 
ture and evolution of such disks see the review of 



Park fc Ostrikerl (|200lh . a nd for the ba s ic phy sics 
of accretion disks see e.g. iFrank et al.l (|2002l ). A 



widely used model f or th e AD is the model of 
IShakura fc Sunvaev ( 1973), with the radial scal- 
ing described in detail in Novikov fc Thorne 
hereaf t er NT), which was also u sed by iRauchl 
(ll995l ): IVilkoviskii fc Czerirvl (j2002h . 

For the radial profile of the AD the surface den- 
sity E is given by 

E(R) = E d f^-\ with a = 3/4, (2) 



y 2 , Rd is the cut-off radius of 



the disk and Ed is the surface density at the cut- 
off radius. The value a = 3/4 corresponds to the 
outer disk range of NT. The mass Md of the disk 
is 



f Rd 2n 

M d = 2tt / Y,(R)RdR = E d i?; 

Jo 2 - " 



(3) 



For the numerical integration using the 4*' l order 
Hermite scheme we need a smooth force in space. 
Therefore we introduce a continuous but steep 
outer cutoff by 



EOR) 



R_ 

Rd 



exp 



-A 



R_ 

Rd 



(4) 



In order to retain exactly Eq. [3] for the total disk 
mass (with the integral now up to oo) we choose 



r i 



a 



5/(2- 



(5) 



with the Gamma- function T(x). For s — > oo the 

cutoff is discontinuous at i? cu t = -Rd with /ji 1 ^' — > 
1. We will use s = 4 leading to (3 S = 0.70 for 
a = 3/4. In this case the surface density at i? d is 
E(i? d ) = 0.49 E d . 

An inspection of the scaling relations in NT 
shows that in the case of SMBHs self-gravity of 
the AD for the vertical structure becomes impor- 
tant for radii larger than ~ lOOrs- Since we cannot 
resolve the innermost part of the AD we simplify 
the model of the vertical disk structure. We adopt 
a self-gravitating isothermal profile given by 



p s {R, z) 



g(g) 



( ) • (6) 



In the NT model the (half-)thickness h z increases 
with distance according oc i? 9 / 8 . Since this is al- 
tered by self-gravitation and since we cannot re- 
solve the vertical structure of such a thin disc close 
to the inner boundary, we decided to adopt a con- 
stant thickness /i z , taking an appropriate value at 
some intermediate radius of the AD. To simplify 
the equations we define the dimensionless value 
h = h z /Rd, which is also constant. 

The effective sound speed c s (which may be 
dominated by the turbulent motion in a clumpy 
gas) is given by (assuming a vertically isothermal 
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model) 

c 2 s (R) = 



8irGp g (R,0)hl = 4V27rGS/i z 



(2 



a)- — h 

fJ-bh 



Rd 



v. 



(7) 

1 (R), 



where v C irc(-R) is approximated by the Kepler ro- 
tation of the SMBH only. In a thin disk the radial 
pressure support can be neglected. Eq. [7] shows 
that the rotation curve u c i rc is highly supersonic 
("^circ/cs ~ 100 at the outer boundary i?d for the 
parameters used in our simulations: /ibh = 0.1, 
/id = 0.01, h — 10~ 3 ). The Mach number is in- 
creasing inwards. 

The sound speed decreases with increasing ra- 
dius and the stability of the AD is a function of 
radius. The Toomre Q stability parameter is given 
by 



Q 



/ kc s y 



2 - a 



2 jJ-bh 



7T /id 



R 



(8) 

with epicyclic frequency k, which shows that 
heavy and thin ADs are unstable near the outer 
boundary. For the case of Hd/Pbh = 0.1 and h = 
10~ 3 the AD is formally unstable at R > 0.26R d . 
In the framework of our simplified AD model with 
constant thickness we may ignore this instability, 
since it could easily be avoided by adopting an 
increasing thickness with distance R. With eqs. 
[5] and [3] the density distribution of the AD with 
constant thickness (Eq. [5]) is given by 



p g (R,z) 



Ma 



2ttV^ h R\ 
R 

Rd 



R 

Rd 



(9) 



cxp 



-A 



exp 



2h 2 R 2 d 



2.2. Dissipation of stellar kinetic energy 

A detailed theory of the dissipative force work- 
ing on the stars crossing the AD depends on the 
details of differential rotation, density profile, tur- 
bulent motion in the disk and possible resonance 
effects. Stars interact with the AD typically 
many times supersonically before they are finally 
trapped in it. Hence, we restrict our investiga- 
tion to supersonic motion only, neglecting the 
details of the last few passages before trapping. 
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Fig. 1. — Drag force |Fd| in units of [7rr^p g Cg] 
for different sound speeds c s in the AD (full red 
lines). We have chosen Qd = 5 and In A = 13. 
The dashed green line shows the approximation 
given in Eq. [10] for the first case. F'geo (dotted 
blue line) is used in the simulations. 



In this case we can use for the dissipative force 
the geometrical cross section F goo = Qd^^JPgKoi 
(determined by the effective area Qd 7rr * of the 
bow shock with stellar radius r* and Qd ~ 5) en- 
hanced by dynamical friction, which is the gravi- 
tational pull by the over-dens ity in the Mach cone 
due to gravitational focusing (|Ostrikerlll999h . The 
dissipative force Fd can be written in the form 
(|Spurzem et al.ll2004l ) 



-ftrlp s \V Ie i\V ie i 



for Vrd > c s . 



In A 



(10) 



Here v csc is the escape velocity at the surface of the 
star (t; eS c=620km/s for the Sun) and the relative 
velocity V rc i = V* — Vd is the velocity of the star in 
the frame co-rotating with the disc. A corresponds 
to the length of the Mach cone in units of the star 
radius r* leading to a Coulomb logarithm In A ~ 
10-20. 

Fig. Q] shows the dissipative force |Fd| for a 
range of values of c s /v csc (full red lines). The 
dashed green line shows the approximation given 
in Eq. [10] for c s /v csc = 0.2. The contribution to 
the dissipative force by Fg CO is shown by the dot- 
ted blue line. From the strong dominance of the 
dynamical friction for velocities V rc \ smaller than 
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v csc we expect that stars are quickly decelerated to 
V Te i < c s and onto co-rotating circular orbits which 
then move slowly to the center with a radial decay 
speed comparable to c s . The dissipative force is 
anti-parallel to the relative velocity and can also 
be accelerating with respect to the rest frame of 
the CSC. 

For a measure of the efficiency of the dissipa- 
tive force we first introduce a dissipative time scale 
^diss*, which describes the energy change pi sd ' of 
a single star due to the dissipative force Fd arising 
from interaction with the AD at the outer part of 
the disk. Our ansatz is: 



tdh 



(11) 



with the stellar mass to* and the 3D stellar ve- 
locity dispersion cr*; the dissipative time scale as 
defined above depends strongly on radius. At the 
outer edge Rd of the AD we have for example 

E* = Q&Krlpgal, p g = Ed/iJd, t dyn = Rd/cr*, 
and one gets 



^diss*(Pd 



0, 



Qd-Pd 



V ' ^dyn 
^d 



(12) 



In the second form of the above equation we 
have defined the surface density of stars = 
TO*/(7rr^), which provides useful insights. Note 
that this dissipative time scale related to a single 
star is a strongly increasing function with radius, 
so the longest dissipation time can be found at the 
outer edge of the disk, in our case at Rd- 

We consider now the quantities £fc and Pd- 
The former implicitly is defined by 2E± := 
£fcTO*erJ(i?d) = Gm+Md/Rd- The latter, Pd, ac- 
counts for the number of disk passages that a star 
will need before its full kinetic energy is dissi- 
pated, and it also includes an average over the 
orbital parameters of the stars, which gives: 



tfvrh 



(13) 



The average in the equation above is taken over 
the disk crossing events of all stars with the proper 
efficiency function g(e,i,R/p) which depends on 
the orbital eccentricity e, the inclination with re- 
spect to the AD i, and the focal parameter p (the 
properties of g(e,i, R/p) will be discussed in de- 
tail in a follow-up paper). For simplicity we ne- 
glect here the contribution by dynamical friction. 



It can easily be included by replacing Qd with the 
velocity dependent factor (Qd + (v csc /V Te \) 4 \nA) 
in the definition of Pd • 

In energy space the dissipation process leads 
to a stationary flow of stars to highly negative 
values; at some point stars will be absorbed by 
the central black hole, for example by tidal dis- 
ruption at the tidal radius r t . As discussed in 
Vilkoviskii fc Czernvl (|2002l ) the dissipation pro- 
cess will bring stars down into the plane of the 
disk, but two-body relaxation will reheat them 
in the vertical direction. Both effects drive the 
stars inwards to the black hole, where they are fi- 
nally tidally disrupted and accreted. Here we use 
a simple parametrized model; stars are absorbed 
onto the central black hole at some accretion ra- 
dius r acc , for which we demand Rd ^> r acc 3> r t , 
but otherwise treat it as a free parameter, varied 
in our simulations later on. 

For scaling purposes with particle number N 
in the simulations we use the standard half-mass 



relaxation time of ISpitzerl (|1987[ ): 



0.147V 



ln(0.47V) 



where rhaif is the half-mass radius of the CSC. The 
half- mass relaxation time t TX is given in Table Q] for 
a series of galactic nuclei covering the observed 
range of SMBH masses. 

The local relaxation time t rc due to two-body 
encounters is given by ( Binnev fc Tremainelll987l 
Eq. (8-71)) 



tr 



0.34cr 3 



G 2 m s p s In A 
18Gyr ( a 



In A Vlkms 1 



lMr. 



(15) 
/ 1M qP c- 3 

V Ps 



with 1-D velocity dispersion a, mean particle mass 
to s , and density p s of the CSC. This equation can 
be applied at each radius or to the central region of 
a system by using mean values inside some radius r 
and is also called the 'core relaxation time'. Inside 
the gravitational influence radius of the black 
hole (which is in our model of the same order as 
the disk outer radius Rd) we assume a\ cx r -1 and 
the standard Bahcall-Wolf density cusp solution 
for th e stellar density profi le in the CSC is p cx 
r' 7 / A dBahcall fc WolJll976l) . Hence we get for the 
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local relaxation time (cf. Eg. [T6|) i rc oc r 1 / 4 . For 
the galactic nuclei in Table Q] the core relaxation 
time at the influence radius i r c(?"h) is a factor of 
1.5 to 2 shorter than t lx . This shows that 2-body 
relaxation in the central part of galactic nuclei is 
well represented by our confined model with a CSC 
mass only 10 times larger than the SMBH mass. 

For comparison with measurements in our sim- 
ulations, as described later, it is more practical to 
define a global dissipation time scale tdiss and a 
global dimensionless dissipative time scale rj in the 
following way: 



tdi: 



^diss / ^rx 

E k _ M cl al(r &cc ) 



E( sd ) 



2E( sd ) 



(16) 
GM hh M cl 



Here E k , with 2E k = M cl al(r acc ) = GM cl M hh /r acc 
and £ , ( sd ) are the total kinetic energy and total 
energy dissipation rate of all stars at the accre- 
tion radius. An accretion time scale for the black 
hole growth is defined by t acc = Mbh/Afbh, or in 
dimensionless form v = £ a ccArx- We remove stars 
at r acc ; so in a stationary state the energy dissipa- 
tion rate required at r acc to sustain the black hole 
mass accretion rate is ij( sd ) = AfbhGMbh/(2r acc ). 
With equation [TTl we get 



^acc ^acc Ad\ 

V = = 77 = T) 

trx ^diss Ad] 



bh 



bh 



(17) 

Hence the black hole mass accretion rate does not 
depend on the choice of r acc , as will be verified in 
our numerical simulations later. 

3. The Numerical Model 

For the high resolution direct iV-body simula- 
tions of the CSC we used the specially developed 
0GRAPE (= Parallel Hermite Integration with 
GRAPE) code. The program uses the 4-th order 
Hermite integration scheme for the particles with 
hierarchical individual block time-steps, together 
with the parallel usage of GRAPE6 (or nowadays 
the GPU - Graphics Processing Unit) cards 
for the hardware calculation of the acceleration a 
and the first time derivative a of the acceleration. 
The code itself and also the special G RAPE hard- 
ware is described in more detail in Harfst et al 



(|2007l ) . For all the new calculations we use the dif- 



ferent NVIDIA CUDA/GPU hardware which em- 
ulate the standard GRAPE library calltQ. 

Here we mention briefly the most important 
features added to the Hermite scheme of the N- 
body 0GRAPE code in order to model the ram 
pressure dissipative force and the star accretion 
to the SMBH: 

• A dissipative force routine, where we calcu- 
late the acceleration caused by the interac- 
tion with the gaseous disk (Eq. |20| and 
its first time derivative ad- 

• We reduce the time-steps of stars when they 
come close to the disk plane. Otherwise 
stars with big individual time step may miss 
the disk and would not feel the effect of the 
disk at all. 

• A simple accretion scheme of stars onto the 
SMBH, where the stellar mass is added to 
the central black hole if they get inside a 
certain accretion radius. The accretion ra- 
dius is used as a free scaling parameter; re- 
sults for different accretion radii can be used 
to scale to real parameters of galactic nuclei 
(see l4.2j) . Thi s algorithm has been described 



and u sed in iFiestas et al.l (|2011l k iLi et al 
(l2012h . 



• In order to control integration error we add 
a careful bookkeeping of energy changes 
caused by removing stars in process of ac- 
cretion and by the dissipative force of the 
stars-gas interaction. In all our runs the 
total energy error does not exceed 10~ 4 . 

In the simulations we use standard A-body 
units given by 



G = M c \ = 4 1 Rot. I = 1, 



(18) 



where E to t is the total energy of the initial Plurn- 
mer sp here for the CSC withou t a perturbing 
SMBH (|Heggie & Mathieul Il986h . Note that 
in this scaling an increase of the particle number 
takes place at constant total mass; hence the stel- 
lar mass cx 1/N decreases with particle num- 
ber, and it also decreases with respect to e.g. the 
disk mass and black hole mass. 



1 ftp -.//ftp. mao . kiev . ua/pub/user s/ber czik/STARDISK/ 
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Table 1: Some examples for the physical scaling of galactic nuclei. 



Object 


A/bh [M Q ] 


cr c [km/s] 


n> [pc] 


TV 


<dyn [Myr] 


*rx [Gyr] 


Qtot 




Otot 


(8*0 


imax [L© 


M 87 


6.6 x 10 9 


312 


291 


6.6 x 10 1U 


1.5 


6 x 10 b 


2.1 


X 


io- 


-a 


5.7 


X 


10- 


-3 


1 x 10 s 


NGC 3115 


9.6 x 10 8 


230 


78 


9.6 x 10 9 


0.58 


3.4 x 10 4 


4.2 


X 


10" 


-9 


1.8 


X 


io- 


-3 


2.8 x 10 7 


NGC 4291 


3.2 x 10 8 


242 


24 


3.2 x 10 9 


0.16 


3400 


1.5 


X 


ur 


-8 


2.4 


X 


io- 


-3 


9.4 x 10 7 


M 31 


1.5 x 10 8 


160 


25 


1.5 x 10 9 


0.26 


2690 


6.2 


X 


io- 


-9 


4.7 


X 


io- 


-4 


5 x 10 7 


NGC 4486A 


1.3 x 10 7 


111 


4.5 


1.3 x 10 8 


0.067 


68.8 


1.7 


X 


io- 


-8 


1.2 


X 


10- 


-4 


1.9 x 10 8 


MW 


4 x 10 6 


HQ 


1.4 


4 x 10 7 


0.021 


7.2 


5.2 


X 


io- 


-8 


1.2 


X 


io- 


-4 


5 x 10 8 


M 32 


3 x 10 6 


75 


2.3 


3 x 10 7 


0.050 


12.9 


1.5 


X 


io- 


-8 


2.8 


X 


io- 


-5 


2 x 10 8 



Giiltekin et al 



Gillessen et al 



Note. — For the CSC we adopt /ij,h = 0.1, solar-type stars with m* = IMq, r+ = 2.3 X 10 8 pc, Vesc = 620 km/s, = 
and T half = 3rh, and = 5 for the bow s hock size. Most obse rvational data (Cols. 2 and 3) are taken from 
i 2009 1). improved values for M 87 are from iMurphv et al.l | |2011| ). and M bh for the Milky Way (MW) is from 
i 2009). The influence radius (Col. 4) is derived from = GM^/cr^, with the core velocity dispersion a c of the CSC; N is 
the number of stars in the CSC derived from Afbh/^bhi the dynamical and half-mass relaxation timescales (Cols. 6 and 7) 
are derived by Eq. 1141 the physical and scaled cross sections Qtot and Qtot(8fc) (Cols. 8 and 9) by Eqs. 1191 and |21I and the 
maximum luminosity by accretion (Col. 10) assuming L max = O-lM^hC 2 . 



Our model consists of three components: the 
CSC, the SMBH, and the AD. The CSC with ini- 
tial mass M c \ is realized by N particles of mass 
m*. The initial density distribution is a Plummcr 
model, which is only in dynamical equilibrium if 
the gravity of the central SMBH is ignored. The 
system adjusts quickly in a few dynamical time 
scales. The mutual gravitational interactions and 
the dissipative force on the stars in the AD are 
calculated fully by the 0GRAPE code. 

The SMBH with initial mass -Mbh = Mbh-Mci 
is modeled by an analytic Kepler potential fixed 
at the origin. We allow accretion of stars, which 
are effectively captured by the inner part of the 
AD or scattered by two-body relaxation into the 
loss cone, to grow the mass of the SMBH. Physi- 
cally the accretion radius r acc is given by the tidal 
radius r t , where the stars are disrupted, or the 
Schwarzschild radius rs. Both radii are orders of 
magnitude below our numerical resolution. There- 
fore we need to analyze the scaling of the accre- 
tion with decreasing accretion radius. The orbits 
of the stars accreted by the interaction with the 
AD are circularized before being accreted. In con- 
trast, stars accreted by two-body scattering into 
the loss cone are predominantly accreted on hy- 
perbolic orbits. In order to simulate the effect of 
the different eccentricity distribution, we apply a 
second accretion criterion based on the velocity 
of the stars. We define an accretion criterion to 
merge stars with the SMBH using two criteria: 1) 
distance r < r acc and 2) V? < k acc v^ irc . Stars 
well inside the influence radius of the SMBH arc 



moving on (local) Kepler orbits, where the veloc- 
ity is given by = v^- 1IC (2 — r / a) with semi- major 
axis a. In the limit of large /c acc all stars reaching 
'"ace would be accreted. For k acc — 1 stars inside 
r acc are accreted, if a < r acc , i.e. all stars with 
energy below —GM c \/r acc (neglecting the poten- 
tial energy of the CSC) are accreted in one orbital 
time. In all simulations we use fc acc = 1. (Note 
for completeness that in some runs with AD we 
accreted all stars with r < r acc /10 independent of 
their energy. We have tested that the accretion 
rates are not changed significantly by this addi- 
tional accretion criterion.) 

The properties of the AD are fixed by the nor- 
malized mass fid with analytic density distribu- 
tion according to Eq. [6] with a — 3/4, s = 4, and 
constant thickness h — 10~ 3 . We use a Kepler 
rotation of the AD in the potential of the SMBH 
neglecting the gravity of the CSC and pressure 
gradients in the AD. We set the outer cutoff ini- 
tially at i?d = n,. Using a Kepler rotation curve 
underestimates the real circular speed at the outer 
boundary by a factor of 1.4, but this has no signif- 
icant influence on the dynamics of the stars. The 
reason is that the gas density and thus the friction 
force is very small in the outer regions of the AD. 
The part of the AD with significant dissipation of 
energy of crossing stars is deep inside the influence 
radius of the SMBH. Here the approximation of 
Kepler rotation is very good. We have chosen the 
large cutoff radius only in order to avoid another 
free parameter. 

It is helpful for calibration of our models with 



S 



respect to real systems in galactic nuclei to define 
an effective dissipative parameter 



Qtot — Qc 



Nnr'1 



(19) 



Note that Qtot enhances Qd by a factor, which 
describes the dimensionless total dissipative cross 
section of N stars, normalized to the disk area. 
With this definition we can rewrite the original 
dissipative force in Eq. [10] as acceleration 



ad = Fd/m* = -Qt 



(20) 



The local parameters like r*, m*, and Qd, whose 
scaling in terms of A-body units is not straight- 
forward, are transformed in this way to the global 
quantities such as Rd, M c \, and Qtot- 

Due to numerical limitations, direct A-body 
simulations are performed here with particle num- 
bers smaller than that in a real galactic nucleus. 
This leads to the relaxation time being a function 
of the particle number, N, whereas it is a fixed 
quantity for a given nucleus. To correct for this, 
the simulations are run for a multiple of the re- 
laxation time and are then rescaled to the correct 
time units to compare to a physical galactic nu- 
cleus, for this reason, all physical processes, in- 
cluding the star-disc drag, need to be rescaled as 
well. In our models all quantities except Qtot are 
invariant when changing the particle number in 
the simulation. Therefore for calibration to real 
systems we just have to compute the value of Qtot 
and adjust it correspondingly in our model system. 
But this step alone is not sufficient; our model sys- 
tem has a much shorter two-body relaxation time 
than in reality. We want to study numerical sys- 
tems at smaller particle number, which have all 
the same ratio 7/ of dissipation to relaxation time 
scale. Therefore we need to choose Qtot such that 
the same value for 77 is achieved: 

QUN 2 ) = t f^l Qtot (N 1 ) (21) 

Our model aims to discover numerically the sec- 
ular evolution of the coupled star gas system in 
galactic nuclei. Therefore we are interested in 
the interplay between dissipative processes (star- 
gas drag) and two-body r elaxation processes - 
the semi-analytic model of IVilkoviskii fc Czernv 



(|2002h defines the setup in which we want to put 
our numerical models. Since we are not yet able to 
simulate star-by-star, two different effects have to 
be taken into account for this: (i) each simulation 
particle effectively represents many stars, not a 
single star. Therefore we cannot study the individ- 
ual star-gas drag effects here, only the collective 
one; this is shown by the definition of Qtot in Eq. 
19\ which shows how for one simulation particle 
an effective cross section is realized, which corre- 
sponds to the effect of the much larger stellar sys- 
tem we have in mind, (ii) In our system with lower 
particle number the two-body relaxation time has 
a "wrong" ratio to other time scales. Therefore, 
we rescale Qtot in such a way, that the relaxation 
time and the dissipation time have the correct ra- 
tio for the larger system, which is done in the Eq. 
2T1 Eq. [21] shows, for example, that for a sys- 
tem with 10 s particles, to be simulated by only 
A = 10 4 , we need to increase Q to t artificially by 
a factor of ~ 5, 000. In Table Q] the real values 
of Qtot and simulated values for A = 8000 are 
both given for the selected galaxies. It should be 
stressed that we do NOT intend to model star-gas 
interactions in detail. Rather we apply a standard 
model of it, scale it up in a collective way and cal- 
ibrate it with respect to the relaxation time. Our 
main goal is to check the increase of star accretion 
rates due to the presence of a central AD. 

In our starting configuration we added the cen- 
tral SMBH to the Plummer distribution of the 
CSC. In order to start the analysis of the system 
in dynamical equilibrium we let the CSC evolve 
for ~ 5tdyn with the SMBH. After the initializing 
phase the cusp around the SMBH shows a density 
slope 7 2.5 which still differs from a BW cusp 
(top panel of Fig. [2]). The difference of the set-up 
Plummer profile and the cumulative mass profile 
at t = demonstrates the effect of virialization 
in the potential of the SMBH. Then we determine 
Rd such that at t = the enclosed CSC mass at 
R = Rd equals the SMBH mass. 

We performed a series of simulations combining 
different particle numbers N = 4fc, 8k, 16k and 
different accretion radii r acc = 0.04, 0.02, 0.01, 0.005 
(see Table E} until t = 10t rx . The identifiers of 
the simulations in Table [5] are coding the particle 
number and accretion radius. For each parame- 
ter combination a run without dissipative force 
is done for comparison. For the different A we 
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Table 2: Model parameters of all runs. 
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-4 
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10 




1 n— 2 
10 


04K-02 


4A- 


4 
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10 


-4 


2 
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04K-04 
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10 


-4 


4 
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io- 2 


08K-005 


8A: 


4 


55 x 


10 


-4 


5 


X 10" 


3 


5.47 x 10- 


08K-01 


8 k 


4 


55 x 


10 


-4 




lO" 2 




5.47 x 10- 


08K-02 


8A: 


4 


55 x 


10 


-4 


2 


x 10- 


2 


5.47 x 10- 


08K-04 


8A: 


4 


55 x 


10 


-4 


4 


x 10~ 


2 


5.47 x 10- 


16K-005 


16k 


4 


55 x 


10 


-4 


5 


x 10- 


3 


2.97 x 10- 


16K-01 


16k 


4 


55 x 


10 


-4 




lO" 2 




2.97 x 10- 


16K-02 


16k 


4 


55 x 


10 


-4 


2 


x 10- 


2 


2.97 x 10- 


16K-04 


16k 


4 


55 x 


10 


-4 


4 


x 10- 


2 


2.97 x 10- 



Note. — Column 1 gives the identification label of the model, columns 2-4 are number of particles, smoothing length and 
accretion radius. Column 5 gives the total cross section, Qtot, scaled according to Eq. 1211 Common parameters for all models 
are fi hh = 0.1, = 0.01, h = IO" 3 , fc aC c = 1. 



applied the scaling according to equation [5TJ For 
the simulations with dissipative force our choice 
of Qtot corresponds to the case of M 87 (compare 
Tables Q] and [5]) . From the last column in Table 
[1] it is obvious that for lower mass black holes the 
value of Qtot hi our simulations is unrealistically 
large. 

For the calculation of the CSC dynamics we 
used the parallel GRAPE systems built at the As- 
tronomisches Rechen-Institut in Heidelberg, and 
at the Fesenkov Astrophysical Institute in Almaty. 
The code has recently been ported to large clus- 
ters with GPU hardware (in Beijing, China and 
Heidelberg, Germany, see acknowledgments) and 
results from these facilities have partly been used 
for this project. 

4. Results 

In our starting configuration the CSC is in 
dynamical equilibrium including the gravitational 
potential of the SMBH and has a steep cusp in 
the density profile. The cumulative mass profiles 
in the top panel of Fig. [2]) shows that after one re- 
laxation time the BW cusp is in place and remains 
stable over the full simulation. The lower panel of 
Fig.[2]shows the Lagrange radii of the CSC for 1, 5, 
25, 50, 90, 99 percent enclosed mass with respect 



to the total mass M c \(t) for the model 16K-005 
(see Table [2] for model parameters). Additionally 
the position of the innermost particle (lp) shows 
that stars crossing r acc on highly eccentric orbits 
are quickly circularized and accreted. The influ- 
ence radius of the SMBH is increasing during 
the simulation by an order of magnitude due to 
accretion and the size of the Bahcall-Wolf cusp is 
growing accordingly. 

4.1. Scaling with particle number 

In simulations without AD, stars reaching r acc 
are immediately accreted onto the SMBH. The 
corresponding loss-cone in phase space has a max- 
imum angular momentum for each energy defin- 
ing the opening angle of the loss cone in phase 
space. The regime of an empty loss cone de- 
pends on tm/tdvn and the width of the loss cone 



2 GRACE: http : //www . ari . uni-heidelberg . de/grace 



(jAmaro-Seoane fc Spurzemll200ll ). In the regime 
with a high probability of a star to be scattered 
in or out of the loss-cone, the loss cone is full. 
In the empty loss-cone regime all stars scattered 
by two-body encounters into the loss cone are ac- 
creted and the accretion rate scales with t rx . In the 
full loss cone regime the accretion rate depends on 
tdyn- In a given system the empty part of the loss 
cone increases with increasing particle number 
leading to a A dependence of the accretion rate. 
We used an additional energy criterion accreting 
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Fig. 2.— Evolution of the CSC for model 16K- 
005. Top panel: Cumulative mass profiles in units 
of the initial CSC mass at different times with the 
Plummer profile and the power law of a Bahcall- 
Wolf (BW) cusp added for comparison. Bottom 
panel: Lagrange radii enclosing a fixed mass frac- 
tion with respect to M c \(t). The position of the 
inner-most particle (lp) and the influence radius 
rh are also shown. 



stars with semi-major axis a < r acc , which cuts 
the loss-cone in energy space. 

In Fig. |3] the lower set of thin lines show the 
growing mass of the SMBH for different particle 
numbers for the cases without dissipation due to 
the AD (top panel: r acc = 0.04 and bottom panel: 
''ace = 0.005). The accretion rate per relaxation 
time is seen to slowly increase with particle num- 
ber. The upper set of thick lines in Fig. [3] show 
that accretion is significantly larger when the dis- 
sipative force of the AD is included. The accretion 
rate is found to be independent of particle num- 
ber N. In Fig. |4] the ratio of the accretion rates 
with and without dissipative force is quantified 
for the 8fc simulations at different accretion radii. 
After a few relaxation times an equilibrium with 
an enhancement factor of ~ 4 is established for all 
accretion radii. 

4.2. The accretion radius 

The physical accretion radius r acc is much 
smaller than the numerical resolution. Therefore 
the scaling of the accretion rate with r acc is very 
important. Without the AD the standard loss 
cone becomes thinner with decreasing r acc , lead- 
ing to a decreasing acc retion rate in the limit of an 
empt y loss cone (e.g. lAmaro-Seoane fc Spurzem 
20011 ). On the other hand the Bahcall-Wolf cusp 
is characterized by a constant mass and energy 
flow to the SMBH and high binding energy, re- 
spectively As a consequence the accretion rate 
based on an energy criterion should be indepen- 
dent of r acc . In Fig. [5] the independence of the 
accretion rate in terms of the accretion timescale 
v (Eq. [TTjl on r acc is shown for all N = 16k runs 
with dissipative force. For t > 1 i rx the accretion 
rate is also independent of the particle number N. 

With the choosen parameters for the AD and 
the CSC the growth timescale of the SMBH by 
accretion of stars is of the order of the half-mass 
relaxation time t rx of the CSC, which is long 
compared to the Eddington accretion timescale 
of m 50Myr (see Table 1). Therefore our model 
applies to quiesent galactic nuclei. The accretion 
rate can be converted to a maximum luminosity by 
adopting L max = O.lMbhC 2 . This upper limit is for 
all SMBH masses in the range of 2 ... 50 x 10 7 L Q 
(last column of Table 1). 
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Fig. 4. — Ratio of the accretion rates with and 
without dissipative force for 7V=8k and different 

^*acc • 
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Fig. 3. — SMBH growth with (upper thick lines) 
and without (lower thin lines) dissipative force for 
different accretion radii r acc = 0.04 (top panel) 
and r acc = 0.005 (bottom panel). The particle 
number ranges from iV=4k...l6k. Q% % scales ac- 
cording to Eq. [21] with N. 




Fig. 5. — Normalized accretion timescale v (Eq. 
[T"7|) for N — 16k with different r acc (same line- 
styles as in Fig. HJ). 
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4.3. Energy dissipation 



Our measurement of the global normalized dis- 
sipative time scale r\ in the simulations uses the 
definition as given in Eq. 1171 

In a stationary state the time scale of transport 
of stars through the AD towards the black hole is 
dominated by the outer edge of the disk, where 
the dissipation time scale is longest. Therefore, a 
measurement of the total dissipated energy in our 
system should not depend on the choice of r acc , 
where we actually remove the stars and add their 
mass to the central SMBH, as it has been shown 
by Eq. 1171 The normalized dissipation timescale r\ 
is a measure of the dissipated energy. It is deter- 
mined by the numerical evaluation of -Ekin(i) and 
E( sd \t), the latter via smoothing of the cumula- 
tive function £A sd )(t) to derive the slope. After 
an initial adaption phase rj is approximately con- 
stant (Fig. IS]) showing the constant energy flow in 
a stationary BW cusp. As seen in Fig. [6l the nor- 
malized dissipation timescale is independent of N 
and depends only weakly on r acc . The long-term 
trend of increasing r\ is due to the evolution of 
the SMBH and the CSC by the accretion of stars. 
The SMBH mass is increasing and the mass of the 
CSC is decreasing leading to an increasing /ibh and 
virial factor fk, and a decreasing total cross section 
Qtot and efficiency with time. By comparing 
the dissipation timescale 77 shown in Fig. |6] and the 
accretion timescale v shown in Fig. [5] we observe 
that the relation derived in Eq. [17] is approxi- 
mately fulfilled (taking into account the growth of 
the black hole and the decreasing CSC mass for 
Mbh)- 

5. Discussion and Conclusions 
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In galactic nuclei super-massive black holes co- 
exist with a dense stellar cluster; galaxy merg- 
ers and the quasar phenomenon indicate that at 
least for some time there should be also large 
amounts of interstellar gas present in the nu- 
clear regions around the black holes. In this pa- 
per we have examined the interaction and co- 
evolution of a dense star cluster surrounding a 
star-accreting super-massive black hole with an as- 
sumed central gaseous disk. Interactions of such 
disks with the surrounding dense star clusters have 
been prop osed as source of gas supply to the cen- 
tral disk (jMiralda-Escude fc Kollmeierl 120051 ) , as 



Fig. 6. — Normalized dissipation timescale r\ (Eq. 
[TTj) for N — 16k and different r acc (top panel) and 
for r acc = 0.005 and different N (bottom panel). 
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agents to enhance the tida l star accretion rate 
(jVilkoviskii fc Czernvl 120021 ). and to cause feed- 



back on the orbital parameters of stars (jRauch 
1995T ) . including a modificati on of sources of grav- 
itational waves (lRauchlll999h . 

Our model is the first self-consistent long-term 
simulation of a dense star cluster, surrounding a 
star-accreting super-massive black hole, and sub- 
jecting the stars to the dissipative forces from a 
resolved central gaseous disk. We resolve effects 
of two-body relaxation, dissipation of stellar ki- 
netic energy in the disk and star accretion onto 
the central black hole in a numerical study based 
on direct high-accuracy simulations. Since star- 
by-star modeling of a galactic nucleus down to 
the realistic tidal radius is not yet possible, de- 
spite of the modern GPU hardware used for sim- 
ulations, a careful scaling analysis is presented as 
a function of the particle number in the simula- 
tions and of an assumed star accretion radius, to 
allow conclusions for the real astrophysical situa- 
tion in galactic nuclei. But our model still has a 
number of serious drawbacks. Firstly, it assumes 
a stationary accretion disk, so energetic feedback 
to the disk structure by star-disk interactions is 
neglected, secondly, the physics of star-gas inter- 
actions is modeled approximately. Finally, there is 
no distinction between properties of different stars 
(main sequence, giants, remnants), but rather a 
single stellar species with solar radius is assumed. 
In that sense our study should be considered as a 
pathfinder and exploratory. 

This investigation is a d irect continuation of 
a sem i-analytic study by IVilkoviskii fc Czernv 
((2002J), extending it by a more detailed and nu- 
merical study of the stellar dynamical evolution 
of the central stellar cluster. Our paper uses a 
numerical approach based on direct ./V-body sim- 
ulations, including particle-particle forces as well 
as a dissipative force in the disk. Here we resolve 
the dissipation of stellar kinetic energy along the 
stellar paths in the vertically extended accretion 
disk. Particle numbers in our simulations and 
an adopted star accretion radius (onto the cen- 
tral super-massive black hole) are used as free 
parameters in our model, while for other impor- 
tant parameters of the problem fiducial values are 
assumed. These are e.g. the initial super-massive 
black hole mass (10 percent of the initial central 
stellar cluster mass), the gaseous accretion disk 



mass (10 percent of the initial black hole mass), 
and the outer radius of the accretion disk (set 
equal to the black hole gravitational influence ra- 
dius). Finally, all star particles are equal in the 
simulation, and their total effective cross section 
is used as a parameter to obtain the physically 
correct ratio of dissipation to relaxation time. 

We show that the accretion rate of stars onto 
the super-massive black hole is strongly enhanced 
by the dissipative force of the accretion disk (a 
factor of four with our parameters). The accre- 
tion process is determined by an equilibrium of 
diffusion by two-body encounters and energy loss 
by the dissipative force. Consistently there is also 
an energy deposition in the central accretion disk; 
we find that most stars accreted or trapped in the 
disk are quickly accreted also to the black hole, 
because there is no stable co-rotation of the stars 
with the disk. Our results are robust with re- 
spect to variation of particle number and adopted 
accretion radius, therefore they should hold un- 
der realistic conditions in galactic nuclei. Our 
star accretion rate does not depend strongly on 
the ad opted star accretion radius; i t supports the 
idea of IVilkoviskii fc Czernvl (|2002l) , that there is 
a stationary flow of stars within the disk towards 
the central black hole, which is determined by 
an equilibrium between dissipation and relaxation 
time. In spite of the enhanced number of stars ac- 
creted through the disk, we still find that there is 
a Bahcall-Wolf central density cusp present in the 
system, which is not significantly perturbed. 

Central densities in star clusters near super- 
massive black holes can reach 10 8 M Q pc -3 or 
more. At such high stellar densities direct, dis- 
ruptive stellar collisions may produce gas in the 
gravitational potential well of the black hole. 
The gas production rate could be larger than 
the one obtained from t idal disruption of star s 
(ISoitzer fc Saslawi Il966t ISpitzer fc Stond Il967tl : 
see al so Begelman fc ReesT i 19781 ); Frank fc Rees 



1976 ), and numerical models in e.g. iFreitag fc Bend 

2002 ); lFreitag. Giirkan fc Rasiol (|2006l) ; lFreitag. Rasio fc Baumgardt 
20061 ). But there is little doubt that a large frac- 



tion of this gas is finally accreted to the SMBH; 
some fraction of it though may escape. The same 
is true for the gas produced by tidal accretion. Our 
models do not yet resolve the very central regions 
of the star cluster, where stellar collisions may 
occur predominantly. We anyway treat the ac- 
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cretion radius as a free parameter used for scaling 
studies; it is here usually large compared to the as- 
trophysically defined tidal radius, where stars are 
destroyed by tidal forces. Our results should not 
depend strongly on whether the stars inside r acc 
are finally disrupted by tidal forces or destroyed by 
stellar collisions with subsequent accretion of the 
gas onto the SMBH. Less is known about the ef- 
fect of induced stellar collisions due to low relative 
velocities of stars in the disk (due to dissipation). 
We will study the effects of direct stellar collisions 
in future work. 

Direct stellar collisions produce another source 
of gas deep in the gravitational well of the cen- 
tral supermassive black hole. As in the case of 
tidal disruptions of stars, we do not know exactly 
how much mass is accreted to the SMBH, and how 
much is ejected due to magnetic fields (jets) and 
radiation pressure near it (see e.g. one improved 
model bv lKasen &: Ramirez-Ruia 2010f) . Our as- 
sumption to add 100% of the material from tidal 
accretion to the black hole clearly is an upper 
limit, the real growth rate may be lower due to 
some mass loss in the process, also for stellar col- 
lisions. However, even in our case with possibly 
overestimated accretion rates, assuming a typical 
value of 10% mass to energy conversion, the lumi- 
nosity obtained from our mass accretion rates are 
of the order of 10 8 L Q , much smaller than the Ed- 
dington luminosity. Therefore our results are ap- 
plicable to quiescent galactic nuclei, not to quasars 
or AGN in their most active phase, where high 
mass accretion rates, feedback and outflows are 
driven for example by gas inflow due to galaxy 
mergers. 

The response of the central accretion disk to 
the deposition of energy and stars is neglected 
in this paper. We will work on this problem in 
future studies. The black hole growth rate due 
to star accretion will be limited by the condition 
that the accretion disk can find a new equilib- 
rium absorbing the dissipated stellar energy and 
radiating it. In the inner spherically symmet- 
ric part of the system the Eddington limit will 
pose an upper limit to the star accretion rate 
allowed in a stationary state (see discussion by 
Miralda-Escude fc Kollmeierll2005| ). In that sense 
our accretion rates, determined by neglecting feed- 
back on the disk, will be upper limits. It is still 
possible that the process of star trapping to the 



disk and star and gas accretion to the central black 
hole is quasi-periodic and highly non-stationary 
(as suggested by the ubiquitous time variability of 
radiation from active galactic nuclei). 
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